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SOME NOTES ON THE NUMERICAL SOLUTION 



OT SHEAR-LAG AND MATHEMATICALLY RELATED PROBLEMS 



By Paul Kuhn 



SUMMARY 



The analysis of "box "beams with, shear deformation, of 
the flanges can he reduced to the solution of a differ- 
ential equation. The same equation is met in other prob- 
lems of stress analysis. No analytical solutions of this 
equation can be given for practical cases, and numerical 
methods of evaluation must he used. Available methods 
are briefly discussed. Two numerical examples show the 
application of the st ep-hy- step method of integration to 
shear-lag problems. 



INTRODUCTION 



Whan a "box beam is subjected " to "bending moments, the 
stress distribution differs somwehat from that given by 
the ordinary theory of "bending . The reason for these dif- 
ferences, denoted by the term tt shear- lag action," lies in 
the fact that the cover sheet suffers appreciable shear 
deformations, particularly after it buckles into diagonal- 
tension fields. Under certain simplifying assumptions 
discussed in reference 1, the analysis of this problem 
leads to a differential equation of the type 

d y a , K 

— ir - Ki y + e s = o (l) 

dx 



In this equation, x is the distance along the span, 
and y may denote, 'for instance, the bending stress in 
the flange, or the shear stress in the cover. The mean- 
ing of Ki and K a depends on the meaning assigned to 
y. 



When a box beam is subjected to torsion, there will 
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be in general a tendency of the cross sections to warp out 
of their plane. At the support, this warping is more or 
less completely prevented and bending stresses arise (ref- 
erence 2). Under certain simplifying assumptions, the 
analysis of this problem also leads to an equation of type 
(1). 

The same equation is obtained in a number of other 
cases. A quite -well-known case is the beam under combined 
lateral loading and axial tension (reference 3). 

Analytical solutions of equation (l) can be given 
when and K a are constants or simple functions of 

x (reference 3). In many practical problems, however, 
Ki and Z a vary along the span in a manner that is dif- 
ficult to define by simple mathematical functions, a dif=- 
ficulty occurring particularly in box-beam problems. It 
is therefore necessary to use numerical methods for the 
solution of given problems. 



NUMERICAL METHODS OF SOLUTION 



All methods to be discussed in this paper depend on 
the assumption that it is permissible to divide the beam 
into a limited number of bays (5 to 10) and that the co- 
efficients K 1 and Z s may be assumed to be constant 
within each bay. This assumption, and the methods used 
for the solution of the differential equation, make the 
solutions inherently approximate. The state of affairs 
is comparable with the engineering usage of applying the 

Mc 

formula <j = -— to tapered beams. 
I 

Trial-a nd- error m ethod.- In reference 4, a somewhat 
unorthodox trial-and-error method was described for ob- 
taining solutions of shear- lag problems. The method was 
found to be very rapid; under favorable circumstances, two 
cycles of the computation were sufficient to obtain an 
accuracy consistent -with the conditions of the problem. 
The speed with which the analysis can be made, however, 
depends critically on the skill and the experience of the 
analyst. It therefore appears desirable to provide other 
methods of analysis that are less dependent on or entire- 
ly independent of the skill of the analyst. Another dis- 
advantage of the trial-and-error method is that it beoomes 
very cumbersome when applied to cases in which the bound- 
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ary conditions differ from those discussed in reference 4. 

F ixed- -po int method. - A graphical method that requires 
no arbitrary trial assumptions whatever has "been described 
Toy L. Kirste. In reference 5, the^method wasepplied to the 
problem of bending stresses dueto torsion; in reference. 6, 
it was applied to the beam-column problem. Reference 6 be- 
ing easily available, the details of the method heed not be 
given here. 

The outline of the procedure is as follows. From one 
end of the beam, the boundary condition given at this Bnd 
being utilized, a sequence of "fixed points" to the other 
end of the beam is constructed. When the other end has 
been reached, the boundary condition given there is uti- 
lized to find the first point of the desired curve, which 
is then established by working back to the original start- 
ing point with the aid of the curve of fixed points. 

Fo r ' purpo ses of comparison with arithmetical methods, 
the drawing of the curve of fixed points may be counted 
as equivalent to one cycle of computation; the drawing of 
the desired curve, as equivalent to a second cycle. Fur- 
thermore, converting the given data into a graphical fig- 
ure and scaling the final curve to obtain the numerical 
answers should be counted as equivalent to one cycle of 
computation. It may be said, therefore, that the graphi- 
cal method is equivalent to a three-cycle arithmetical 
method. 

The neatness and the straightforwardness of the 
fixed-point method are impressive. /It is doubtful," how- 
ever, whether it will be possible in Borne cases to achieve 
the necessary accuracy by graphical means. There' 'are al so 
practical objections to graphical methods on the grounds 
of ease of checking, filing, and transmitting computations. 
For such reasons, graphical computations are usually em- 
ployed only when they save considerable time as compared 
with arithmetical methods. It is very questionable whether 
the fixed-point method has a very marked advantage over 
the step-by-step integration method to be" next described. 

S.tep.- b y- a t ep integration method.- An orthodox step- 
by- step integration method may be used either in semi- 
graphical or in numerical form. Only the numerical form 
will be described here because it is believed to be pref- 
erable for all-around use. 
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It is impossible to start the final integration at 
either end of the team "because only one "boundary condition 
is given at either end. It is therefore necessary to make 
an arbitrary trial assumption for the unknown boundary con 
dition at one end before starting a first integration. 
When this first integration is completed, the "boundary con 
dition at the far end will not have been met (except by 
accident). A second integration is therefore made, start- 
ing with a different initial assumption for the unkown' 
boundary condition and resulting in a different error of 
closure. If the errors of closure are plotted against the 
initial- assumptions , it is possible to find what initial 
assumption must be made to reduce the error of closure to 
zero . 

The graphical form of this method was described in 
reference 7 as applied to the beam-column problem. One 
quite important point, however, was apparently not noted 
or at least was not specifically mentioned in this paper. 
The curve of error of closure against -initial assumptions 
is a straight line. Consequently, two trials are suffi- 
cient, in principle at least, to determine the desired in- 
itial assumption and it will not be" "irecessary to make more 
than three cycles of computation. The third cycle will 
be the final one, a fact that may be used as a check 
against - "erTo"rs of computations. 

In the appendix, detailed examples are given which 
show the application of the numerical procedure to shear- 
lag problems. The basic theory and the nomenclature is 
taken from reference 1. 

Remarks on the ac c uracy of numerical 'solut ion s.- The 
errors of numerical solutions may be divided into three 
classes: 

1. Errors arising from using a small number of sig- 

nificant figures. 

2. Errors arising from using bays of finite length. 

3. Errors arising from the assumption that the coef- 

ficients are constant in a bay. 

Errors of the first class are important only when the cal- 
culation involves small differences between large figures. 
In most cases, this difficulty can be overcome by using a 
calculating machine; ordinarily a slide rule is sufficient 
ly accurate.. 
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A different method of overcoming the difficulty is 
to rearrange the computation so that no Bmall differences 
are involved. In the case of shear- lag problems, this 
rearrangement can "be accomplished by writing the equa- 
tions not for the actual stresses hut for the corrections 
that must "be applied to the stresses given by the ordi- 
nary "bending theory. (See appendix.) This procedure is 
probably more rational than the first one in several re- 
spects. , 

« 

The magnitude of the errors of the second class can 
be estimated by making comparisons with known analytical 
solutions. Tor simple cases, it was found that five or 
six bays often give sufficient accuracy (about 1 percent). 
One rather extreme case investigated was a beam column 
with a compressive load near the buckling load (L/j = 3). 
In this case, the half-beam had to be divided into 10 
bays. (Only the half-beam was used on account of sym- 
metry. ) 

Errors of the third class, can also be appraised by 
making comparisons with analytical solutions; they ap- 
pear to be of the same order of magnitude as the errors 
of the first two classes. 

The conclusion may be drawn, then, that the necessity 
for making "approximate" numerical solutions of the dif- 
ferential equation (l) is not, in general, the factor 
limiting the practical accuracy of the results. The fac- 
tors limiting this accuracy are the simplifying assump- 
tions, physical and mathematical, that must be made to 
reduce the problem to the solution of a simple differen- 
tial equation. 

-n ... 

Langley Memorial Aeronautical Laboratory, 

National Advisory Committee for Aeronautics, 
langley yield, Va. , March 31, 1939. 
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APPENDIX 



APPLICATION OP THE METHOD OP STEP- BY-STEP INTEGRATION 

TO S HEAP- LAG PROBLEMS 



Th e d i fferential equations of the shear-lag theory. - 
Two static equations and one elastic equation form the 
"basic equations of the shear-lag theory for "box "beams with 
flat covers (reference l). 

<V Aj = !l - T t (i) 

h I 



V a l 



= Tt 



t — — — — (a, - <r T .) 



Eh 



(2) 
(3) 



The symbols used are those of reference 1. The prime de- 
notes differentiation with respect to x. The x origin 
will "be taken at the tip of the beam for all cases. 

The "basic equations can he used to form the follow- 
ing differential equations. 

<T « - K a a_ + K 3 JL. = 0 (4) 
7 *A T 

<T L " - E 3 CT L + K 2 =0 (5) 

hArp 

S A 

t" - K 8 t + K 3 -i-Jt =0 (6) 

ht A T 



where 



a & e Wl 1 
E b ^A p A L 



(~ + — ) (?) 
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For convenience, the "boundary conditions are summarized 
in table I. ....... 

It will "be noted that the equations for Oy and 

are identical, so that cry and cr^ can he obtained by 

solving the same differential equation for different 
"boundary conditions at the root. It will "be found, how- 
ever, that the stresses ohtained in this manner Batisfy 
the equation of static equilibrium with the external mo- 
ment only at the root and at such stations along the span 
where the ratio A^/A^ is the same as at the root. If 

this ratio varies radically along the span, the solution 
will "be reliable only near the root. 

The numerical evaluation of the equations requires 
the computation of the second" derivative. Prom equation 
(4), .for example, ' 




This expression may "be very small over a large part of 
the span, making it necessary to carry a large number of 
significant figures. It is therefore better to write the 
differential equation for the expression in parentheses, 
which is the correction that must be added to the stress 
calculated by the ordinary bending theory to account for 
shear lag, because a = K/hAijj 1b the stress calculated 

by the ordinary theory. 

If this correction is denoted by 



(10) 

where Z has the same meaning as before. This equation 
can be solved by step-by-step integration, but the solu- 
tion is somewhat easier than the Solution of one of the 
equations (4), (5), or (6). 



u = 



M 

hA 



T 



the differential equation becomes 

u" - E 2 u = 0 



8 



U.A.C.A. Technical Mote No. 704 



In analogy with equation (9), the correction for the 
stress cr^ may be defined "by 

hA T 

giving the differential equation 

v« - K a v = 0 (10a) 

The solutions of (10) and (10a) for any given problem 
again fail to satisfy the equation of static equilibrium 
with the external moment except at the root and at those 
stations where the ratio Ap/A^ is the same as at the 

root. The discrepancy was always found to be smaller than 
when solving equations (4) and (5). For general U6e, 
equation (10) therefore appears to be the best choice. 
The stresses t and cr-^ can be found from Btatics after 

u, and therfore , have been found. 

The equations discussed thus far apply to box beams 
with flat covers. For box beams with camber (fig. 2), 
equations corresponding to (4), (5), and (6) could be 
written. As in the case of .beams without camber, it ap- 
pears to be more convenient to use equation (10). The 
constant K in equation (7) is then changed to 



Eb« \ A F A L 
Example 1 

Find the stress (T^, in the axially loaded panel shown 
in figure 3. The following data are given: 

P 1 

=1.00 A- = A_ = - A m L = S.00 

Am F L 2 T 



x' . 1.00 
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(These numerical values were chosen to emphasize method 
rather than arithmetic. They do not imply impractical 
sizes; they imply merely unconventional units for stress 
and length, which need not he specified here.) 

The panel will be divided into five bays. Stresses 
and other functions at any station will be designated by 
a -subscript denoting the values of x/L in tenths; for 
instance, cr 4 designates the value of Op at the station 
x/l =0.4. 

Two trial values must be assumed for ffj 1 . If the 
shear stiffness were infinite, the stress would be con- 
stant, and the corresponding value a* = 0 will be used 
as the first trial value. If the shear stiffness is 
large, but finite, the value of cr at the root will be 
nearly 



o- 10 = — = 1.00 



while at the tip 



cr o = = 2.00 



A- 



If 



Assuming tliat the stress decreases linearly between the 
tip and the root gives 

Act 1.00 

cr' = — = = -0. 500 

• Ax 2.00 



This value will be used for the second trial. 



Table II shows the arrangement of the computation. 
The start of the first cycle of computation is as follows. 
At x/L = 0, the known boundary value cr = 2.000 is writ- 
ten down. The first trial value for cr 1 * is zero; there- 
fore, 

i.l- * t-i ■ ■ ■ - =- 

Acr - = cr' Ax = 0 X 0.400 = 0 

O "~ 2 

and consequently 
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Ca = Ofi! + = 2.000 + 0 = 2.000 
Now, using equation (4) in the form 

V = K n S <J n - K n 2 (f-) 

x n 

giveB 

CTa" a 2.000 - 1.000 = 1.000 

(The subscript n denotes the station number and is added 
here to indicate that, in the general case, K and P/Aj 
vary along the span.) 

Therefore 

Acr x _ 3 ' = cr s " Ax = 1.000 X 0.400 = 0.400 

Adding this increment to the value at the preceding sta- 
tion gives 

ct 3 » = <j 1 i + Ao- 1 _ 3 ' = 0 + 0.400 a 0.400 

Prom here on, the computation repeats itself in principle. 
Prom the slope a' just found, the increment Ac is cal- 
culated. Adding this increment to the stress at the pre- 
ceding station gives the stress at the new station. Using- 
the differential equation gives a" , etc. In the general 
case, K and P/A T will be different for each station; 
the table would therefore include a row for K a <7 and a row 
for K a M/hAj. 

Since an axially loaded panel can be considered as 
the cover of a beam in pure bending, the boundary condi- 
tion that must be satisfied at the root is, according to 
table I, .■ r~" 



oW = 



hA 



P 



In the first two trials, the values of ff 10 T obtained are 

2.905 and 1..047, respectively. (Note that these values 
are obtained by adding l/2 Axcr" to c 9 ' . ) Thoso values 
of a* aro obtained by assuming cr^ = 0 and CTj} = -0.500, 
respoctivoly; by proportion, it is found that oV should 
bo takon as -0.782 to satisfy tho condition ct 10 ' = 0. This 

initial value of a* was used in the third cycle of tho com- 
putation; the computation shows that (7 X0 l =0 as required 
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(with a very small error). Comparison of the stresses 
with those found "by formula (roforonco 1 or 2.) shows that 
the error in the root stress is 0.3 percent. 



Example 2 

As a second exampl e , ' the analysis of a tapered "beam 
with cambered cover will he given. The "beam chosen is 
N.A.C.A. "beam 4; the dimensions of this beam, as well as 
the results of strain-gage tests, are given in reference 
1. The methods used in obtaining the "basic data given 
in table 111, .rows 2 to 6, will "be briefly outlined. 

The compression flanges, the tension flanges, and 
the longitudinals are assumed to he concentrated at their 
respective centroids. The cover sheet is assumed to be 
^ully effective in aiding the longitudinals, and the ad- 
jacent strips of sheet are added to the longitudinals or 
flanges. The web is also. assumed to he fully effective 
in bending and is replaced by concentrated flanges of 
cro s s- sect ional area l/6 A^. 'The area of the " substitute 
longitudinal" is calculated by 

sinh K 3 b 
Hs = A L 

K 3 h 

which is formula (4) of reference 1; the "substitute cam- 
ber" is taken as eg = l/z c. The thickness of the cover 

sheet is 0.0114 inch, and the ratio G/E is taken as 
0.40. With the values given in rows 2 to 6, the values 
of k can be computed by formula (7a). 

If an arbitrary value of u' = 1.00 is assumed for 
the first bay, a step-by-step integration of equation (10) 
can then be made, as shown in rows 8 to 12. 

According to table I, the boundary condition that 
must be fulfilled at the root is 

hAy Aj " 

Now at the root 

M 

5 W = P - - tan a = 0.473 P 
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giving 

0;473 P X 0.463 



u i = 



5.70 X 0.177 X 0.640 



For an applied load P of 250 pounds, this equation be- 
comes 

u« = 85.0 

The value obtained in the trial solution is u' = 247.1, 
so that the final values for u are obtained by multi- 
plying the u values of the trial solution by 

3 5 0 

'. — = 0.344. These final u values are given in row 

247 - 1 Mc 

13. In row 14 are given the values of — cr =• — , that 1b, 

I 

the flange stresses obtained by the ordinary bending the- 
ory, using the geometric properties of the beam given in 
reference 1. Adding the u values to the cr values 
gives the final solution as shown in row 15. Figure 4 
gives a comparison between calculated and experimental 
stresses . 
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TABLE I 
BOUNDARY CONDITIONS 



Tip 




"D - -v Ju / -v 




( x — d 1 


Longitudinal 
"built in 


JJU X U U.U. -U net 

not "built in 


Elastically yielding 
support (fig, 1) 


M 
hA p 


^ = hi; 


it 

1 nA F 


1 hA F Ay o 


a L = 0 


a L " = 0 


a L = 0 




G eM 


T = 0 




T = V & Q 


T = - 

EohAp 


T 1 - 

EDhAj, 


M A L 
u =■ — 


u t = S ff A L 
hAj, A T 


hA F Aj, 


u , = JL y a e 

hAj, Ag, A-g, 


h*\ T Aj, 


V = 0 ' 


. S W 
v ' = -— 


_ M 

hArp 


t' = — V f. 

hArp A^ 



r 
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v TABLE II 

SHEAR- LAG- ANALYSIS OF AXIALLY LOADED PANELS 



x/L 


0 




0.2 




0.4 




0.6 




0.8 




1.0 


Firs 


fc trial solution 


a' 




0 




0.400 




0.864 




1.466 




2.303 


2.905 


Aa- 




0 




.160 




.346 




,586 




.921 




CT 


2.000 




2.000 




2.160 




2.506 




S.092 




4.013 


o-" 






1.000 




1.160 




1.506 




2.092 




3.013 


Ao-' 






.400 




.464 




.602 




.837 




1.205 



Second trial solution 



0" 1 




-0.500 




-0.180 




0.111 


l 


0.420 




0.796 


1.047 


ACT 




-.200 




-.072 




.044 




.168 




.318 




cr 


2.000 




1.800 




1.728 




1.772 




1.940 




2.258' 


cr" 






.800 




.728 




.772 




.940 




1.258 


Act' 






.320 




* 

.291 




.309 




.376 




♦ 503 










Final 


solution 










ct' 




-0.782 




-0.507 




-0.313 




-0.169 




-0.053 


o.poi 


Ac- 




-.313 




-.203 




-.125 




-.068 




-.021 




er 


a. 000 




1.687 




1.484 




1.359 




1.291 




l.£70 


ct" 






.687 




.484 




.359 




.291 




.270 


Act' 






.275 




.194 




.144 




.116 




.108 


Analytical solution 


CT 


3.000 




1.684 




1.482 




1.355 




1.287 




1.2 66 



TABLE III 

SESAE-UG AHA1YSIS OF TAPE8ED BEAM WITH 0A1SERED COVER 



1 




u 












O f 










2 




2.70 




3.30 




3.90 




4.50 




5.10 




5.70 




A— ■ ( an \ yi . 1 


.143 




.150 




.157 




.163 




.170 




.177 


A 
*± 


A — f en ^ vi \ 






?72 


















O 


h t fin ^ 






( 


















6 


Co fin.} 
S ^ ' 






.67 




.78 




.90 




1.01 




1.13 


7 


IT 


.00967! 




.00729 




.00575 




.00472 




.00394 




.00334 


8 


u' 




1.00 




3.63 




13.23 




43.68 




131.1 


[247.1] 


q 


An 








fin q 








rro 








10 


u(lbysq.in.) 


0 




19.0 




87.9 




339.4 




1,169 




3,659 


n 


u» 






.1385 




.505 




1.603 




4.60 




12.23 


12 








2.63 




9.60 




30.45 




87.4 




232.0 


13 








7 




30 




117 




402 




1,260 


14 


cr(HVsq.in.) 


0 




1,795 




3,690 




2,920 




3,040 




3,015 


15 


a F t HVsq.in.) 


0 




1,802 




2,720 




3,037 




3,442 




4,275 



o 



H3 
O 

a 

H» 

O 
P 
H 

!z! 
o 

CD 

!zi 

O 



O 



H 
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Fig. 4 
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Figure 4.- Comparison between calculated and experimental stresses 
in 1T.A.C.A. "beam 4. Data from reference 1. 



